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Abstract 



oo 

f" — ■ Tectonic faults are commonly modelled as Volterra or Somigliana dislocations in an elastic 

l' medium. Various solution methods exist for this problem. However, the methods used in prac- 

^— ^ tice are often limiting, motivated by reasons of computational efRciency rather than geophysical 

rf\ accuracy. A typical geophysical application involves inverse problems for which many different 

,-H fault configurations need to be examined, each adding to the computational load. In practice, 

IL" this precludes conventional finite-element methods, which suffer a large computational over- 

. I^ head on account of geometric changes. This paper presents a new non-conforming finite-element 

S^ method based on weak imposition of the displacement discontinuity. The weak imposition of 

;_( the discontinuity enables the application of approximation spaces that are independent of the 

Cu dislocation geometry, thus enabling optimal reuse of computational components. Such reuse 

of computational components renders finite-element modeling a viable option for inverse prob- 
lems in geophysical applications. A detailed analysis of the approximation properties of the 
new formulation is provided. The analysis is supported by numerical experiments in 2D and 3D. 

Keywords: Volterra dislocation. Finite Element Method, weak imposition, linear elasticity, 
tectonophysics. 
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1 Introduction 

The world is perpetually reminded of the fact that seismic hazard is still beyond reach of prediction 
— as it was most recently by the disaster that struck Japan. The difficulty is not just to predict 
the exact moment of failure, which, as argued by some |9], might never reach a level of practicality. 
It is also the nature of the risk, and the extent to which stress is accumulating, that turns out to 
be surprisingly difficult to constrain. The 2011 Tohoku-Oki earthquake demonstrated a great lack 
of understanding of ongoing tectonics [12] . Arguably, a better understanding could have reduced 
the secondary effects if such information had led to more apt measures and regulations. 

The main reason for this poor state of information can be traced to the absence of direct 
measurements. The primary quantities of interest, being the magnitude and orientation of the 
stress tensor in the earth's crust, can be obtained only through tedious, expensive, point-wise 
measurements. A viable broad scale method to directly measure the global stress field does not exist. 
For this reason information is obtained mostly from secondary observables, earthquakes themselves 
being an important source. Earthquakes represent significant, near instantaneous changes in the 
global stress field. By accurately determining the location of the segment of the fault that collapsed, 
one can progressively update the stress field and evolve it in time. This way the tectonic evolution is 
monitored, and hazardous areas can be identified as regions where stress accumulates. For successful 
tracking of stress development, however, it is essential to understand the tectonic mechanism behind 
every earthquake. This includes the location and geometry of the section of the fault that collapsed, 
and the direction and magnitude of fault slip. It is increasingly popular to base such analyses on 
local co-seismic surface displacements. This type of information has become available since the 
nineties with the advent of space borne interferometric SAR measurements of the earth's surface, 
and with the widespread availability of GPS measurements 19! . Analysis of this data has in recent 
years seen rapid adoption and is now routinely performed for all major earthquakes. 

A mechanical model is required to connect observations to physics. Most commonly (if not 
exclusively) used is an elastic dislocation model, based on the assumption that on short time scales, 
nonlinear (plastic) effects arc negligible. The model embeds a displacement discontinuity of given 
location and magnitude in an elastic medium, causing the entire medium to deform under the 
locked-in stress. Many different solution methods have been developed for this particular problem, 
based on analytical solutions or numerical approximations; see for instance 21j for an overview of 
the most prominent methods. However, methods founded on analytical solutions generally dictate 
severe model simplifications, such as elastic homogeneity or generic geometries, which restricts their 
validity. The computational complexity of methods based on numerical approximation, on the other 
hand, is typically prohibitive in practical applications. Because in practice the surface displacements 
are given, and the dislocation parameters are the unknowns, the computational setting is always 
that of an inverse problem. A typical inversion requires several thousands of evaluations of the 
forward model, and therefore computational efficiency is a key requirement. Moreover, the forward 
problems in the inversion process are essentially identical, except for the fault geometry. Reuse 
of computational components, such as approximate factors of the system matrix, is imperative for 
efficiency of the inversion. Current numerical methods for seismic problems do not offer such reuse 
options. 

Finite-element methods provide a class of numerical techniques that are particularly versatile in 
terms of modeling capabilities in geophysics. Finite-element methods allow for elastic heterogeneity, 
anisotropy, and topography; all things that can not well be accounted for with currently used 
analytical and semi-analytical methods. In geophysical practice, finite-element methods are however 



often rejected for reasons of computational cost. The high computational cost can be retraced to 
the condition that the geometry of the fault coincides with element edges, which is a requirement 
engendered by the strong enforcement of the dislocation; see |l4j. Consequently, the mesh geometry 
depends on the fault, which in turn implies that mesh- dependent components such as the stiffness 
matrix and approximate factorizations of that matrix cannot be reused for different fault geometries 
and must be recomputed whenever the geometry of the fault changes. The recomputation of these 
components in each step of a nonlinear inversion process leads to a prohibitive overall computational 
complexity. 

To overcome the complications of standard finite-element techniques in nonlinear inversion pro- 
cesses in tectonophysics, this paper introduces the Weakly- enforced Slip Method (WSM), a new 
numerical method in which displacement discontinuities are weakly imposed. The WSM formula- 
tion is similar to Nitsche's variational principle for enforcing Dirichlet boundary conditions [T5] . 
The weak imposition of the discontinuity in WSM decouples the finite element mesh from the ge- 
ometry of the fault, which renders the stiffness matrix and derived objects such as approximate 
factors independent of the fault and enables reuse of these objects. Therefore, even though the 
computational work required for a single realisation of the fault geometry is comparable to that of 
standard FEM, reuse of components makes WSM significantly more efficient when many different 
fault geometries are considered. This makes finite-element computations based on WSM a viable 
option for nonlinear inverse problems. 

A characteristic feature of WSM is that it employs standard continuous finite-element approxi- 
mation spaces, as opposed to the conventional FEM split-node approach |14j which introduces actual 
discontinuities in the approximation space. Instead, WSM approximations feature a 'smeared out' 
jump with sharply localized gradients. We will establish that the error in the WSM approxima- 
tion converges only as 0{h^^^) in the X^-norm as the mesh width h tends to zero and that the 
error diverges as 0{h^^'^) in the energy norm. In addition, however, we will show that the WSM 
approximation displays optimal local convergence in the energy norm, i.e., optimal convergence 
rates are obtained on any subdomain excluding a neighborhood of the dislocation. The numerical 
experiments convey that WSM also displays optimal local convergence in the L^-norm. 

The remainder of this manuscript is organized as follows. Section[2]presents strong and weak for- 
mulations of Volterra's dislocation problem, and derives the corresponding lift-based finite-element 
formulation. Section[3]introduces the Weakly-enforced Slip Method based on two formal derivations, 
viz., by collapsing the support of the lift onto the fault and by application of Nitsche's variational 
principle. In Section l4l we examine the approximation properties of the WSM formulation. Sec- 
tion [5] verifies and illustrates the approximation properties on the basis of numerical experiments 
for several 2D and 3D test cases. In addition, to illustrate the generality of WSM, in the numerical 
experiments we consider several test cases that violate the conditions underlying the error estimates 
in Section |4j such as discontinuous slip distributions and rupturing dislocations. Section [6] presents 
concluding remarks. 

2 Problem formulation 

In this section we define Volterra's dislocation problem. We will postulate the strong formulation 
in |2.1[ followed by a derivation of the weak formulation in |2.2| The latter will serve as a basis for 



the derivation of the Finite Element approximations, for which we lay foundations in Section 2.3 




Figure 1: Definition of the domain fi, normal vector v, slip vector 6, fault V and dislocation ><. 

2.1 The strong form 

We start by defining the geometric setup. We consider an open bounded domain Q, C M^ {N = 2, 3) 
with Lipschitz boundary dVt. An (A^— l)-dimensional Lipschitz manifold F, referred as the fault, 
divides the domain in two disjoint open subdomains rj"*" and VL~ , such that i7 = int(ri+ U r2~). 
We equip F with a unit normal vector v : T —^ E^ directed into the subdomain D,~ . The fault 
supports a slip distribution 6 : F — > M^, corresponding to a dislocation. The fault is referred to as a 
non-rupturing fault if the dislocation x — supp(6) is compact in F, and a rupturing fault otherwise. 
Let us note that in tectonics, rupturing faults correspond to intersections of the slip plane with the 
surface of the earth. Figure [T] illustrates this setup for N — 2. It is to be noted that b need not be 
tangential to the fault. If b has a non-vanishing normal component then the fault is opening, such 
as may be caused by an intruding material. 

The displacement field generated by the dislocation is represented by u : f7 \ F ^- M^. For 
convenience, we restrict our considerations to linear elasticity. We denote by the map u n- e(ii) the 
strain tensor corresponding to the displacement field u, according to 

e(w):=i[Vw+(Vuf], (1) 

under the assumption that u is differentiable on fJ \ F. The constitutive behavior corresponds to 
Hooke's law: 

a(w):=A:e(u), (2) 

i.e., (Jij{u) = Aijkiekiiu), where we adhere to the convention on summation on repeated indices. 
The tensors a and A are referred to as the stress tensor and the elasticity tensor, respectively. The 
elasticity tensor is subject to the usual symmetries Aiji^i = Aijn- = A^nj. Moreover, we assume 
that it is bounded and satisfies a strong positivity condition, i.e., there exist positive constants 
CA > and qa > such that: 

for all tensors e. The elasticity tensor is in principle allowed to vary over the domain il, subject to 
the above conditions. Auxiliary smoothness conditions on A will be introduced later. 

To facilitate the formulation, we denote by cr„(u) := cr(u) ■ n the traction on a boundary corre- 
sponding to the displacement field u. Moreover, we introduce the jump operator [[•]] and average 
operator {•} according to: 

lv]\ -.r 3 X ^ v+{x) - v^{x), (4a) 

{v} : r 3 X ^ ^[v+ {x) + V- (x)] , (4b) 



where w"*" and v^ represent the traces of v from within Q^ and D,^ , respectively. Given a partition 
of the boundary dfl into T) ^ % and M such that 2? n A/" = 0, we define the Voherra dislocation 
problem as follows: 
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Equation (5a I is the usual equilibrium condition, which applies everywhere in f2 except on the 



manifold V. Equations (5b) and (5c) respectively express that the displacements at the boundaries 
of il+ and 51^ differ by the slip vector b, and that the tractions at the boundaries of fi"*" and 
rj^ are in static equilibrium, i.e., equal and opposite. It is to be remarked that this condition 
corresponds to a standard linear approximation in the small-slip limit, as the traction equilibrium 



at the fault occurs in fact in the deformed configuration. The boundary conditions ( 5d ) and ( 5e I 
correspond to Dirichlet and Neumann boundary conditions, representing a prescribed displacement 
and a prescribed traction, respectively. 

To facilitate the presentation, we note that the solution to the Volterra dislocation problem (l5| 
can be separated into a discontinuous part uq : Vi\T — > M^ with homogeneous data and a continuous 
part ui : 17 — >■ M with inhomogeneous data: 



(6) 



The sum ug + ""i satisfies (|5]). Therefore, the inhomogeneous data /, g, /i in ([5]) can be treated 
separately in a standard continuous elasticity problem, and without loss of generality we can restrict 
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our consideration to homogeneous data. We retain / := to identify the right member of (5a). 

For rupturing faults, some compatibility conditions arise with respect to the boundary condi- 
tions. In particular, an intersection of the dislocation with the boundary of the domain is only 
admissible at the Neumann boundary M . Otherwise, an inadmissible incompatibility between the 
jump condition [[wq]] = & and the homogeneous Dirichlet condition mq = ensues. 



2.2 The weak form 

To derive the weak formulation of ([5]), we note that for any piecewise smooth function v from fi \ F 
into M^, we have the identities: 

V ■ divCT(u) = IX) V ■ an{u) — / cr(u) : Vw 
o\r Jan- Jn- 

V ■ (Tn{u) — / cr{u) : Vu 
aa+ Jn+ 



V ■ cr„(w) — / o'(u) : Vv 
on Jn\r 

+ lM-WAu)} + {v}-[{aAun (7) 

The first identity resuhs from integration by parts. The second identity follows from a rearrange- 
ment of the boundary terms and 

'"^ ■ O-n (w) + «" • 0-r7(w) = V+ ■ cr+(u) - V^ ■ CT^ (u) 

r Jr 

{v+-v-)^{a+{u) + a-{u)) 

+ Jj{v++v-){a+{u)-a-{u)) (8) 

In the weak formulation, the admissible displacement fields will be insufficiently regular to ensure 
the existence of the tractions cr„(-). Hence, the terms involving these tractions in ([t]) must be 
eliminated by means of the boundary conditions and auxiliary conditions on v. The traction term 



on the Neumann boundary can be eliminated by means of ( 5e ). To remove the traction term on dVl, 
we stipulate that v vanishes on V. The traction average in the final term of ^ is eliminated by 
requiring that v be continuous. The traction jump in the final term is deleted by means of the 



traction-continuity condition ( 5c ) . 

Summarizing, wc find that a solution u of ([5]) satisfies 

ar{u,v) = l{v) (9) 

for all sufficiently smooth functions u : J7 ^- M.^ that vanish on V, where 

ar{u,v) — / a{u):Vv (10a) 

Jn\r 



Kv) ^ v-f (10b) 

Jn 

Note that v is assumed to be smooth on il and, in particular, that it is continuous across the fault 

r. 

To furnish a functional setting for the weak formulation of (l5| based on ^, we denote by 
H'^ (ri) the usual Sobolev space of square- integrable functions from ft into M.^ with square- integrable 
distributional derivatives of order < fc, equipped with the inner product 



and the corresponding norm || • \\k,n and semi-norm | • \k,n- For the square-integrable functions and 
the corresponding norm and inner-product, we introduce the condensed notation L'^{Q) := ff°(r2), 
II • ||n := II • ||o,n and (•, •)n := (•, •)o,n- We denote by i?J ^(O) the subspace of H^{fl) of functions 
that vanish on I? C dfl. To accommodate the discontinuity corresponding to the dislocation, we 
introduce the hft operator £(.), which assigns to any suitable slip b : T —>■ M.^ a function ii, in 
JEfj!)^(ri \ r) such that [[£(,]] = b. A precise specification of conditions on the slip distribution is 
given in Section l4] The weak formulation of ^ based on Q writes 



Weak formulation: given the lift £& e Hq ^(51 \ F), find u E l^ + H\ j,{^) such that 

ar{u,v) = l{v) \/veHlj,lSl). (11) 



The bilinear form ar : H^[n\V) x H^{n\T) -> M and linear form I : H^{n\T) ^ M in ( 11 ), are the 



extensions (by continuity) of the corresponding forms in (10). The treatment of the dislocation by 



means of a lift operator in ( 11 ) is analogous to the treatment of inhomogeneous Dirichlet boundary 



conditions in weak formulations; see, e.g., [7| p. 113 ]. Let us note that in the weak formulation ( 11 ), 
we have identified {u € H^{U, \ F) : [[u]] = 0} with ff^(fi). This identification is unambiguous, on 
account of a one-to-one correspondence between the functions in these spaces. 



To analyze the weak formulation ( 11 ), and to prepare the presentation of the Weakly-enforced- 
Slip method in Section pj we note that (11) is to be interpreted in the following manner: find 



u :— u + ti, with u G JfJ p(ri) such that 

a{u,v)^l{v)-aT{ib,v) yv & Hlj,{n) . (12) 

where the bihnear form a : H^{Q) x if ^(51) -^ M, 

a{u,v) = / (t{u):\/v, (13) 

Jn 

corresponds to the restriction of ar(-, •) to H^{rt) x H^{n). Indeed, for all pairs {u, v) E H^{il) x 
H^{n,), the function cr(u) : Vv is Lebesgue integrable on fl, and because the manifold F corresponds 
to a set of A^-Lebesgue measure zero, the integrals of (t{u) : Vv on J7 \ F and on il coincide. It 



is important to note that the restriction of the bilinear form ar(-, •) to H^{fl) x H^{il) in (13) is 



independent of the fault F. The function m in (11) is referred to as the continuous complement of 



the solution u with respect to the jump lift £i, and, indeed, it resides in Hq •p(ri). 

For the assumed linear-elastic behavior according to (IT]) and (l2|, it follows straightforwardly 
that 

\ar{u, v)\<CA |w|i,n\r|i^li,n\r (14a) 

|ar(w, u)| >CA |w|i,o\r (14b) 

for all u,v € H (J7 \ F), where ca and ca denote the continuity and strong-positivity constants 
of the elasticity tensor, respectively, and | • |i_a\r represents the usual H -seminorm. By virtue of 
Poincare's inequality (see, e.g., [2 Theorem 5.3.5]) there exists a bounded positive constant Cp 
such that 

\\u\\i,n\r < Cp\u\,^n\r Vu e if ^^^(r! \ F) (15) 

7 



Hence, the H -norm and H -semi norm are equivalent on Hq •p(ri). Equations ( 14) and ( 15 1 imply 
that the bilinear forms ar(-,-) and, accordingly, a(-, •), are continuous and coercive. Moreover, it 
is easily verified that the linear form l{-) — ar{ib, •) : H (S7) — >■ K in the right member of (12) is 



continuous. Problem (12) therefore complies with the conditions of the Lax-Milgram lemma (see, 



for instance, [H Theorem 2.7.7]) and, hence, it is well posed. 



It is interesting to note that (12) allows for a physical interpretation that is very close to 
Volterra's classical construction for dislocations, popularly known as the 'Volterra knife': to make 
a cut in the material, displace the two sides and hold them while welding the seam, and finally 
release the sides so the material assumes its state of self-stressed equilibrium. The initial cut and 
displacement is represented by the lift i, which is not in equilibrium and is hence maintained by 
an external load. The addition of u represents the transition to a state of equilibrium, by removing 
the external load but leaving the displacement intact. 

2.3 Finite element approximation 

Galerkin finite-element approximation methods for Volterra's dislocation problem ( |11[ ) are generally 
based on a restriction of the weak formulation to a suitable finite dimensional subspace. The general 
structure of finite-element methods can be found in many textbooks, for instance, [Tl] [531 [71 13] . We 
present here the main concepts and definitions for the ensuing exposition. 

The approximation spaces in finite-element methods are generally subordinate to a mesh Th, 
viz., a cover of the domain by non-overlapping element domains k G fl. The subscript h > 
indicates the dependence of the mesh on a resolution parameter, for instance, the diameter of the 
largest element in the mesh. In general, we impose some auxiliary conditions on the mesh, such 
as shape-regularity of the elements and conditions on the connectivity between elements; see, for 
instance, [71[S] for further details. A finite-element approximation space V^ C Hq p(fi) subordinate 
to Th can then be defined, for instance, as the subspace of vector-valued continuous element-wise 
polynomials of degree < p which vanish on V: 

VI = {vh e C"(n,M^) : {vhUU e PP for all k e Th,i ^ 1,2, . . .,N,Vh\v = 0} (16) 

with PP the 7V-variate polynomials of degree p. Below, our interest is generally restricted to the 
/i-dependence of the approximation space and, accordingly, we will suppress p. The finite-element 



approximation of (11 ) based on an approximation space Vh writes: find u/i :— Uh+it with Uh € Vh 
subject to 

a{uh,Vh) = l{vh) ~ar{£b,Vh) Vw/^ G F^ . (17) 



We refer the right-most term in the right member of (17) as the lift term. 

Approximation properties of the Finite Element Method are generally investigated on the basis 
of a sequence of meshes Tn '■= {Th)h£'H, parametrized by a decreasing sequence of mesh parameters 
T-L — {/ii, /i2, . . . } with as only accumulation point. A sequence of meshes is called quasi uniform 
if there exist positive constants C and C, independent of h, such that Ch < diam(K) < Ch for all 
K G Th and all h G H. Standard interpolation theory in Sobolev spaces (see, for instance, [71 [2]) 



conveys that a sequence of approximation spaces V-^ of the form ( 16 ) based on quasi-uniform meshes 
possesses the following approximation property: there exists a positive constant "^^ independent 

tWe use V to denote a generic positive constant, of which the value and connotation may change from one instance 
to the next, even within a single chain of expressions. 



of h such that for ah h £ H, all k > Q and both m E {0, 1}, it holds that 

inf \\v-Vh\U,n<'^h'+^-"'\v\i+,^n Vv e H'^+^n) D H] ^{n) (18) 

Vh&Vh 



with I = min{p, A: + 1}. The estimate (18) imparts that for all sufficiently smooth v e ffQp(51), 



the 11 • llm^Q-norm of the best approximation in Vh in that norm decays as 0(ft-^+^^™) as /i — )■ 0. 

The lift £b in ( |17[ ) is in principle arbitrary. However, the use of an arbitrary lift carries severe 
algorithmic complexity, as one has to explicitly construct the lift and evaluate integrals involving 
products of (gradients of) the lift and finite-element shape functions. Moreover, the evaluation of 
these integrals by a suitable numerical integration scheme generally leads to a high computational 
complexity, because the fault is allowed to intersect elements, and there are no efficient quadrature 
schemes to integrate the discontinuous function that arises. Therefore, in practice, it is convenient 
to integrate the lift in the finite-element setting. Provided that that the fault coincides with element 
edges, a lift Ith is then constructed in the broken approximation space: 

Vh = {ue C"(n\ r) : u\n± £ (Vh)b±} (19) 

It is to be noted that the slip does not generally reside in [[V/i]] and, accordingly, b is to be 
replaced by a suitable interpolant bh- Moreover, if the fault does not coincide with element edges, 
then it is to be replaced by an approximation subject to this condition. The aforementioned 
approach corresponds to the split-node method by Melosh [T3], where the adjective 'split' refers to 
the discontinuities between the elements in £7"*" and 0~ contiguous to F. The split-node approach is 
analogous to the standard treatment of Dirichlet boundary conditions; see, for instance, [71 §3.2.2]. 
The split-node approach bypasses the aforementioned complications of an arbitrary-lift approach 
and the evaluation of the lift term comes essentially free of charge as part of the regular stiffness 
matrix. A fundamental disadvantage of the split-node approach, however, is that it requires that 
the fault coincides with element edges, which connects the fault geometry to the geometry and, 
generally, the topology of the mesh. As a result, computational primitives such as the stiffness 
matrix and preconditioners for the stiffness matrix, which are contingent on the mesh, cannot 
be reused for analyses of alternative fault geometries. This is a prohibitive restriction if many 
dislocation geometries have to be considered, for instance, in inverse problems. 

3 The Weakly-enforced Slip Method 



In section [273] we substantiated that the treatment of the lift term in the split-node approach, which 
is the natural counterpart of the standard treatment of Dirichlet boundary conditions in finite- 
element approximations, is unsuitable if many fault geometries have to be analysed, on account of 
the inherent dependence of the finite-element mesh on the fault geometry. 

In this section we propose a new and fundamentally different treatment of the lift term that 
retains mesh independence: the Weakly- enforced Slip Method (WSM). Below, we present two dif- 
ferent formal derivations of the WSM formulation. The derivation in Section 13.11 relies on a limit 
procedure. In Section |3.2[ we derive the WSM formulation on the basis of Nitsche's variational 
principle for enforcing Dirichlet-type boundary conditions ^15) . 



3.1 Collapsing the lift 



Our aim is to derive a tractable finite-element approximation of Volterra's dislocation problem ( 11 ) 



in which the finite- element space and the bilinear form and, accordingly, the stiffness matrix are 
independent of the fault geometry. 



In principle, the lift-based Galerkin formulation (17) already exhibits the appropriate form 



However, as elaborated in Section 2.3 the corresponding finite-element formulation is intractable 
for general lift operators. One can infer, however, that the complications engendered by a general 
lift operator can be avoided by collapsing the support of the lift on the fault. The integration of a 
discontinuous function in H. then reduces to the integration of a smooth function on the dislocation. 
Numerical evaluation of integrals on the dislocation is feasible given a parametrization of the fault. 
Moreover, the intricate explicit construction of a lift is obviated, and only the slip distribution 
itself is required, which is presented as part of the problem specification. A further advantage of 
collapsing the lift is that of localization: instead of having to evaluate the lift term for all shape 
functions of which the support intersects with the support of the lift, only the shape functions of 
which the support intersects with the dislocation have to be considered. 

To derive the lift term corresponding to a collapsed lift, we consider a symmetric lift £| as 
illustrated in Fig. [2j with compact support in an e-neighborhood of the fault: 



{x en :dist(j;,r) < e}. 



By virtue of the compact support of l^ in T^ , the lift term of (17) evaluates to 



ariilv) = 



a{v):V£l= / b-{a,iv)}- / il ■ div a{v) . 
r=\r Jr Jr^\r 



(20) 



(21) 



The first identity follows from the symmetry of the bilinear form in ( 10a ) . The second identity 



results from integration-by-parts and [[^|]] = b and {£1} — 0. Let us note that the second identity 
is formal in the sense that it requires more regularity of v than is actually provided by H^{n,). We 
shall momentarily ignore this aspect, but it manifests itself in the analysis of the approximation 
properties of the WSM formulation in Section l4] Without loss of generality, we can assume il to 



be bounded independent of e. The second term in the final expression in (21 ) then vanishes if T^ 
collapses on F. Therefore, formally passing to the limit in (21), we obtain 



ar{ei,v) 



>+o. 



b-{a.{v)} 



(22) 



According to ( 22 ) , the lift term reduces to an integral on F in the limit of collapsing the support of 
the lift onto the fault. The WSM formulation corresponds to replacing the lift term ar{£b, •) in the 



right member of (17) by the limit functional according to (22): 
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Figure 2: Schematic representation of an e-local lift supp(£j) G F^, which is added to a continuous w to 
form the discontinuous displacement field u. 



The nomenclature Weakly- enforced Slip Method serves to indicate that in (23) the slip discon 



tinuity is weakly enforced in the right-hand side only, and does not appear in the approximation 
space. 



It is to be noted that although the WSM formulation ( 23 1 is derived from the lift-based formula- 



tion (17) by collapsing the lift, in contrast to (17) we do not add a lift to the continuous complement 
Uh- Because Vh C C°(0,]R^), WSM thus yields a continuous approximation to the discontinuous 
solution of the Volterra dislocation problem (11). This implies that the approximation near the 



dislocation will inevitably be inaccurate. We will however show in Section [4] that away from the 
dislocation, the error in the WSM approximation converges optimally under mesh refinement. 



3.2 Alternative derivation via Nitsche's variational principle 

To further elucidate the WSM formulation, we present in this section an alternative derivation 
of ( 23 ) based on Nitsche 's Variationsprinzip (15j . Nitsche presented in |15j a variational principle 



for weakly imposing Dirichlet boundary conditions in finite-element approximations of elliptic prob- 
lems, i.e., without incorporating such essential boundary conditions in the approximation space. 
Nitsche's variational principle can be extended to the Volterra dislocation problem (111 to weakly 



impose the slip discontinuity. To specify this extension, we consider a suitable broken space V{h) 
which encapsulates the broken approximation space Vh and contains the solution u to the Volterra 



dislocation problem (11). We define the quadratic functional J : V{h) 



J{w) = -ar(w, w)- [[w]] ■ {cr,.(w)} -I- 



w 



(24) 



for some suitable constant -0 > 0, generally dependent on h. Let Vh denote either the broken 
approximation space Vh or the continuous approximation space Vh and consider the approximation 
Uh^ u according to: 

Uh'-= aig'vc£ J{u-Vh) (25) 

VheVh 



Equation (25) implies that Uh satisfies the Kuhn- Tucker optimality condition J'{u — Uh){vh) = 
for all Vh e Vh^ where v k> J'{w){v) denotes the Frechet derivative of J at w. For J according 
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to (24), the optimality condition implies that Uh S V^/i satisfies: 

ar{uh,Vh)- [[uh]] ■ {cTuivh)} - [[vh]] ■ {cru{uh)} + i^ I [[uh]] ■ [[vfi] 



^ar{u,Vh)- I M ■ {a,{vh)} - J [[vh]] ■ {<J,{u)} + i^ / [[«]]•[[«]]/, (26) 



= l{vh)~ j b ■ {a,{vh)} + ^J J b-[[v]]h yvheVh. 

The final identity follows by invoking integration-by-parts on ar(M, w/j.), a rearrangement of terms 
and the strong formulation of Volterra's dislocation problem in (l5|. 

If the broken approximation space Vh is inserted for V^, the optimality condition (26 1 can be 
reinterpreted as a symmetric-interior-penalty (SIP) discontinuous-Galerkin-type formulation; see, 
for instance, Sec. 4.2]. In contrast to standard discontinuous Galerkin formulations, the slip 



terms in the right-hand side, i.e., the terms containing b in the ultimate expression in (26 1, enforce 



the jump discontinuity at the fault. Convergence results for this formulation can be established 



in a similar manner as in [15!. For suitable stabilization parameters ip, the functional J in (24) is 



equivalent to || • ||i,r2\r ^^^d (25) implies quasi-optimal convergence of m/j. 



If the continuous approximation space Vh is inserted for Vh, the terms containing [[u/i]] and 



[[tift]] vanish, and we obtain the WSM formulation (23). Hence, WSM can indeed be interpreted as 
an extension of Nitsche's variational principle to the Volterra dislocation problem with continuous 
approximation spaces. Furthermore, in view of Vh 3 Vh, WSM can also be regarded as a SIP 
discontinuous Galerkin formulation, based on a continuous subspace. One can infer that the WSM 
approximation retains the quasi-optimal approximation property in || • ||i_n\r- However, since the 
continuous approximation spaces applied in WSM are not dense in HQ-i~,{fl \ F), the immediate 
significance of this quasi-optimality for the approximation properties of WSM is limited. 

4 Approximation properties of WSM 

An analysis of the approximation properties of the Weakly-enforced Slip Method is non-trivial, 
owing to the fact that in WSM one considers approximations in Jifp p(r2)-conforming subspaces, 
while the solution itself resides in Hq p(ri \ F), and the embedding of Hq •p(^) in Hq p(f2 \ F) 
is non-dense. Essentially, we attempt to approximate a discontinuous function by a continuous 
one and, in doing so, we incur an error that does not vanish under mesh refinement. Standard 
techniques to assess global approximation properties, on all of il \ F, viz., Cea's lemma or the 



Strang lemmas [7l Lems. 2.25-27], therefore provide only partial information; see Section 4.2 

To provide a foundation for analyzing the approximation provided by WSM, we we first recall 
some aspects of traces and tractions in Section [4T| Section [42| investigates the global approximation 
properties of WSM, i.e., on the entire domain. Section [4.3| establishes the local approximation 
properties of WSM, i.e, on the domain excluding a neighborhood of the fault. 

4.1 Traces and tractions 

To enable an analysis of the approximation behavior of WSM, some elementary aspects of trace 
theory are required. For a comprehensive overview, we refer to [18j . To make the theory applica- 
ble to the Volterra dislocation problem, we must impose auxiliary smoothness conditions on the 
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elasticity tensor. In particular, we assume: 

A,kieC\n) (27) 



It is to be noted that ( 27 1 implies that the elasticity tensor is C^ continuous on the domain, including 
the boundary, and across the fault, including the dislocation. The C^ continuity on the subdomains 
0+ and n~ ensures that tractions are well defined. The C^ continuity across the fault is required 
to establish global convergenc e of the WS M approximation in the X^-norm and local convergence 



in the if ^-norm; see Sections 4.2 and 4.3 below. Let us note that in tectonophysics the elasticity 
tensor is generally assumed to be uniformly constant in the domain. 

Let to C K^ denote an arbitrary connected domain with Lipschitz boundary. In particular, 
recalling the partition of H, into the complementary subsets $7^, we envisage ui S {57+, 57"}. We 
denote by 7 the restriction of a function in C^{uj) to the boundary dui. By virtue of the density of 
C (w) in H (w), the operator can be extended to a linear continuous trace operator, denoted by 7, 
from H (w) into L (duj). The image of 7 is denoted by if ' (duj). Considering a subset x C duj, 



we denote by 7^ := (7(-))U the composition of the trace operator and the restriction to xr. The 
image of ^ 
Hq (>f): 



image of j^ restricted to the class of functions if q ouiX^ti'^) that vanish on du} \ >;: is indicated by 



Hl/\>,) = {j^u:ueHlg^^^{u:)}. (28) 

1/2 
The space if q (>f) can be endowed with the norm: 

||A||i/2,:^ :=inf {||u||i^^ : u e Hlg^\^{uj),j^u = X}. (29) 

with the obvious extension to if ^"(9w). There exist continuous right inverses 

of 7 and 7>f. Such a right inverse is called a lifting (or lift) of the trace. It is to be noted that lift 
operators are generally non-unique. 

We denote by ct„ : w 1— > n • 7(cr(u)) the traction of a function in C {uJ) on duj, where n denotes 
the exterior unit normal vector on duj. We define 

Hl,,,iuj) := {u e H\uj) : diva(w) e L^u;)} . (31) 

Applying index notation for transparency, the chain rule yields 

djcrij{u) = {djAijki)<^ki{u) + Aijki{djeki{u)) (32) 

Therefore, the condition div cr(u) G L (u) provides a meaningful condition on u if djAijki G L°°{lu). 
For uj G {ri"'"j_£2^}, this auxiliary condition on the elasticity tensor is satisfied under the standing 



assumption (27 1. The vector space H^^^^^uj) is a Hilbert space under the inner-product associ- 
ated with the norm (|| • Hii^ + ||divcr(-)||^)^/^. The traction (?„ can be extended to a bounded 
linear operator, denoted by cr„, from H\^^^{u}) into H^^''^{duj) :— [if ^"(9(u;)] , the dual space of 

if ^"(9aj). For each u G H]^i^^{uj), the functional crn(u) acts on functions in if ^"(9aj) by means 
of the following duality pairing: 

(a„(M),A)= [ divaiu)-^-\X)+ [ a{u) :Vj-\X) (33) 
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One may note that for functfons in C (cj), Equatfon (33) corresponds to a standard integration- 
by-parts identity. Continuity of the operator cr„ thus defined follows from the sequence of bounds: 



|(cr„(w),A)| < ||divcr(u)|| ||7"^(A)|| +|k(u)|| |7"^(A)L 

< (lldivaHll^ + llaHliy (||7-^(A)||: + |7-^(A)| J 
<(l + ci)'/^(IHlL + ||divaHP^/^||7-^(A)||^^ 



1/2 



(34) 



and the continuity of the lifting of the trace from H.^^'^ {dui) into H^{uj). The dual space [iJ ^ " (t^"^)] 
is a Banach space under the norm 



ll^ll-l/2,cl<. 



sup 

Aeffi/2(aa 



{v,\) 



lAlli 



/2,9a 



The restriction of the traction ((?„(•)) |>c to a subset >f C duj of the boundary can be extended 
to a bounded linear operator an >t from H^-^^^{uj) into H~ ' (k) :— [Hq {x)] ■ The functional 
Cr»,j«r(w) acts on functions in H^ [x) via the duality pairing: 



(a„^^H,A) = / di^a{u)-^-\\) + / a{u) : V7;;^(A) 



(35) 



Continuity of the operator cr„ j^ thus defined follows in a similar manner as in ( 34 ) 

.1/2 



For non-rupturing faults, it holds that h G Hq (F) and the above definitions apply without 
revisions. The slip discontinuity ( 5b ) and the traction discontinuity ( 5c ) are then to be understood in 



the sense of traces and tractions outlined above. However, for rupturing faults, i.e., if the dislocation 
intersects with the boundary of the domain, then h 
We can accommodate h in the space: 



1 /2 

if Q (F), and further consideration is required. 



if i/2(F) := {7rM : u e H^{uj)] 
with u) e {5^^, 51^}. The space ii-'^/^(F) is a Banach space under the norm 

l|A||i/2,r ^ inf {||u||i,c.. : u e ii^(w), 7rM = A} 



(36) 



(37) 



The principal complication pertaining to rupturing faults, is that the corresponding slip vectors 
cannot be lifted into iio,9tj\r('^)i ^s traces of functions in H (w) do not admit the discontinuity 
that would otherwise arise at the intersection of F and dio \ F. Hence, we cannot use (35 1 to define 
an extraision of the restriction of the traction, (CTn(-))|r, to a bounded lineax_operator from H^:^^^{u]) 
into [ii^/^(F)] . However, there exists a continuous right inverse 7^^ : ii^/^(F) -^ H (cj) of the 
operator 7r, for instance. 



7p^(A) = arginf ||w|i,j^ : u e H^{uj),jru = x\ 



(38) 



Let us note that the image of the lift operator 7p corresponds to a harmonic function subject 
to inhomogcneous Dirichlet conditions on F with data A and a homogeneous Neumann condition 
on duj \ F. The lift operator 7^ can be modified to include homogeneous Dirichlet conditions on 
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V C duj \ r in the codoniain, if necessary. The hft operator 7p^ enables us to extend (CT„)|r to a 
continuous hnear operator: 

5„,r : [u e -ffdiv^(w) : \\crn,au\r{u)\\_^/^g^^^ = o| ^ [H^^^{r)]' (39) 

The functional (Jnxi''^) ^'^^^ o^^ functions in H^''^(r) via the duality pairing: 

(a„j(u),A> - / diva(M) -TpHA) + / (riu) : VTr^A) (40) 



Essentially, in (39) and (40 1, we have defined the extension an,r of the restriction of the traction to 
the dislocation, (CT„(-))|r, by restricting the domain of the extended operator to functions for which 
the traction vanishes on duj \ T. This restriction in the definition is consistent with the standing 
assumption that an intersection of the dislocation with the boundary of the domain can only occur 
at Neumann boundaries. 



In the analysis below, we restrict ourselves to non-rupturing faults. The analysis in Section 4.2 
however extends to rupturing faults by replacing the spaces and trace and traction operators for 
non-rupturing faults with those for rupturing faults. 

4.2 Global approximation properties of WSM 

To assess the global approximation properties of WSM, we first construct an upper bound on the 
functional (6, {cr^{-)})r ■ V h — ?> M in the right member of the WSM formulation. It is to be noted 
that, in general, Vh <(- -ffdivcr(^)- Hence, the functional (6, {cri,(-)})r : V y^ — > M does not admit an 



interpretation as a duality pairing according to ( 35 1 . Because 'V h is piecewise polynomial, however. 



an upper bound can be constructed on the basis of inverse and trace inequalities. We refer to [5, for 
a comprehensive treatment of this subject. Inverse and trace inequalities can generally be derived 
under suitable (sufficient) regularity conditions on the finite-element mesh; see [51, Chap. 1]. A 
detailed treatment of the conditions underlying inverse and trace inequalities is beyond the scope 
of this work. Instead, we shall directly assume that a suitable discrete trace inequality holds. To 
formulate the assumption, for a given sequence of partitions Th and for each /i e "H, we denote by 
Sh a dense cover of the fault by means of open intersections of the fault with element interiors or 
with element faces, i.e., 

5,i = {s C r : s = r n K 7^ for some k e %} 

U {s C r : s = int(r n Okq n Oki) ^ for some kq, ki G Th, kq ^ ki} 

See Figure[3]for an illustration. The separate treatment of element boundaries serves to ensure that 
subsets of r that coincide with element faces are separately included in Sh- To each segment s d Sh, 
we associate a pair of contiguous elements {k+, kJ} such that s C k* U dn^ and k^ il^ ^ 0. 
If s C Ok (resp. s C k) for some k then k^ and k^ will be distinct (resp. identical). We assume 
that the following discrete trace inequality holds for all ft. G "H, all s G Sh and all element-wise 
polynomial functions vt of degree at most p: 

{di^m{K)Y^^\\vh\\s<Cr\\vhh, (42) 

for both K G {k+, kJ}, for some Cr > independent of h; cf. [5, Lemma 1.46]. The constant Cr is 
allowed to increase with the polynomial order p. 
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Figure 3: Illustration of the covering of the fault T by segments Sh = {a, b, c, b, e}, and the corresponding 
elements {/ct, ,k7, }. Because b,c,d coincide with element boundaries nf, 7^ nT, for these segments. The 
segments a, e are interior to elements and, accordingly nf, = k7. for these segments. 



Lemma 1 (Continuity of the WSM linear form). Consider a manifold F C il C M.^ , a slip 
distribution b G Hq (T) and a sequence of partitions T-h such that for all h GH, the discrete trace 



inequality [^2] holds for all s £ Sh and all element-wise polynomial functions on Th and 

1/2 



\ThX 



E'^^^'ii^ii' <^ 



seSh 



with hs the harmonic average of the diameters of the elements adjacent to s, 

1 _ 1 1 

hs diam(KJ' ) diani(K7) 

Then for all h&'H, the linear form {b,{(Ju{')})v ■ Vh — > K is continuous and 

|(6,K(t;)})r| < 2-1/2 c^Crikfi/^ |j^|j^^^^ ||^||^^^ 



(43) 



(44) 



(45) 



with ca the continuity constant of the elasticity tensor in Q), Cp the constant in the discrete trace 
inequality (4-2) and M the maximum multiplicity of the multiset {n £ {k+, kJ} : s € 5^}. 



Remark. The maximum multiplicity M in Lemma [T] indicates the maximum number of occurrences 
of any one element in connection to any segment s g 5^ as a member of the set {nfjKj}. For 
instance, in Figure ^ the element k+ = k" = k^ has multiplicity 3. One can infer that M is 
bounded by the maximum number of faces of any clement in the mesh, increased by 2 for interior 
segments. 

Proof. We first separate the integral on F into a sum of contributions from the segments and 
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apply ([3]) and the Cauchy-Schwarz inequality to obtain 



|(&, {cr^(uh)})r| 



2-1 Y. (^^•7i^('TK)) + ;^-7r(^K))), 



seSh 



<2-'cA J2 ll^ll-(l|7r(Vf/^) 



|7r i^'^h) 



seSh 



(46) 



where 7p (•) denotes the trace of (•) from within J7*. Noting that Vh G Vh is clement-wise polyno- 



and (44) 



mial, we deduce from the discrete trace inequality (42), the arithmetic- geometric mean inequality 

M|Vt;,.|L.,) 



|(6,KK)})r| < 2-ic^Cr ^ h-'/'\\b\\,{\\Vvh\\^- 

seSh 

<2-'/'cACrM'/^\\b\\nxhhKn 



1/2 



(47) 



D 



To determine the global approximation properties of the WSM formulation in the H -norm, 
we note that the WSM formulation is inconsistent: The solution of the weak formulation of the 



Volterra dislocation problem (11 ) violates the WSM weak form (23) by (&, {a^{v)})r, for all v e Vh 



The second Strang lemma ^ Lemma 2.25] then provides the following characterization of the global 
approximation properties of WSM: 

Theorem 2 (Global approximation properties of WSM in the H -norm). Assume that the con- 
ditions of Lemma \l\ hold. Let u € iJQp(J7\r) denote the solution to the Volterra dislocation 
problem ITTp and let Uh € Vh denote its WSM approximation according to {23). It holds that 



\u-Uh\\i.Q,\T <{^ + Cpc/ca) inf ||u - Wft|li.n\r 



+ Cpc^ sup 



l(&,K(^h)})r| 
VheVt,\{o} \\vh\\i,n 



(48) 



Proof. We first recall that the bilinear form ar(-, •) is bounded on H\ ^{Q. \ F) x Vh and coercive 
on Vh X Vh] see (14) and (15). For arbitrary Vh £ Vh, we have the chain of inequalities: 



W'^h -f^/i||?,o\r < Cpc^^\ariuh -Vh,Uh - Vh)\ 

= Cpc^^\ariuh -u,Uh- Vh) +ar{u-Vh,Uh - Vh)\ 
< Cpc^^ {\{b, {ai,{uh - Vh)})r\ +ca\u - Vh\i,n\r\uh ~ Vh\i,n\r) 
which leads to 

m{aAwh)})r\ 



(49) 



\uh - Vh\\i,n\r < Cpc/ sup 

wheVh\{o} 



\wh\\i,n 



+ Cpc^ ca\\u ~ Vh\\i^n\i 



The bound ( 48 ) then follows from the triangle inequality. 



(50) 

D 
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For quasi-uniform meshes, LemmafTland Theorem [2] lead to a simple asymptotic characterization 
of the global approximation properties of WSM in the H^-norm. This characterization is detailed 
in the following corollary: 

Corollary 3. Assume that the conditions of Lemma^^hold and that the sequence of partitions Th 
is quasi-uniform with respect to the mesh parameter, i.e., for all h Cz H there exist constants C > 
and C > independent of h such that Ch < diam(K) < Ch for all k CzTh- It then holds that 



\\u-uh\\i,n\r<'rfh-'^^ 
as ft, — > 0. 
Proof. Subject to the quasi-uniformity condition on the sequence of partitions, we have 



(51) 



l^'llr. 



r 



seSh 



1/2 



< 2C-i/2ft-i/2 



(EIIHI^)'^' = 2C:-V2/,-i/2|,;,|,^ (52) 



It then follows from (47) and (48 1 that 

\W~Uh\\i,n\r < (1 

< 



inf |lM-u/,.|li,f2\r 
vheVh 



CpC^ CA 

2-i/2Cpc^icACrMi/2|15|lr,,r 



and the assertion follows in the limit /i — > 0. 



(53) 



n 



The potential divergence of the bound in Corollary [3] does of course not immediately imply that 
the error in the WSM approximation itself increases as h — > 0. However, Theorem |4] below asserts 
that \\u — Uh\\n vanishes as ft — > 0. From this result, the continuity of Uh and the discontinuity of 
u, one can infer that, indeed, \\u — M/i||i.fi\r rnust diverge as ft — > 0. 

Theorem I2I conveys that for each ft > 0, it holds that u — u^ e L'^{fl). Based on the Aubin- 
Nitsche Lemma, we can can therefore construct an estimate of the error in the WSM approximation 
in the L -norm. The estimate is formulated in Theorem W^ below. Some auxiliary conditions on the 
domain fl are required, as specified in the premises of the theorem. 

Theorem 4 (Global approximation properties of WSM in the i^-norm). Assume that the con- 
ditions of Lemmaln hold. In addition, assume that fl is convex or of class C^ . Let u denote the 

solution to (uTp and let Uh G Vh denote the WSM approximation according to (23). Let Th denote 

— 1/2 — 

the intersection of^riVh) with Hq (F). It holds that 



|(u-u/,.,i^)| < ''#'( |w- u/i|io\r + ll^llr^.r) inf \z^ - Vh\i.Q 

\v\\n\\b- Xh\\ -\ 



'if inf 

AhGTh 



(||A/i||i/2,r + i|A/i||rfc,r) inf \z^ - Wh\i,n\r 

Wh&Vh 



(54) 



for arbitrary Lp £ L (Jl) and a corresponding z^ Cz H (fl) D Hq^{CI) such that \\zip\\2.n ^ ^Hv'lln; 
with II • llj-^^r defined by (43). 
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Proof. For arbitrary ip (£ L (il), consider the dual problem (in the sense of distributions): 





—diva{z^) = ip in ri 




z^ = on P 




Cniz^) =0 on TV 


or, equivalently, in weak form: 




z e Hlj,m : 


ar{v,Zip) = {v,f)n 



yveHlj,{n) 



(55a) 
(55b) 
(55c) 



(56) 



By virtue of the smoothness conditions on the elasticity tensor (27), the conditions on the domain. 



and the smoothness of the (homogeneous) boundary data in (55b I and (55c), the dual problem (55) 



possesses an elliptic-regularity property (see, e.g., [HI [J)- The regularity property implies that the 
dual solution resides in H'^{Q) n Hlj,{fl) and satisfies the estimate ||z^||2,a < '^||</'||n- 

Denoting by if, g JEfJp(ri \ T) a suitable hft of b such that [[4]] = b, it holds that u ~ ib & 



HQj,{fl). Moreover, the WSM approximation Uh resides in ijjp(r2). The dual problem (56) 
therefore gives: 



{u - Uh, <f)n = {u~ 4, ip)i-i + (4, ip)i-i - (u,j, Lp)i-i 

^ ariu - 4, z^) + (4, ip)n - ariuh, z^) 

= ar{u, z^) ~ ar{u,i, z^) - (ar(4, z^) - (4, <p)o) 



(57) 



By means of ( 33 ) , the term in parenthesis in the ultimate expression can be identified as the weak 
formulation of (5, {crjy(z^)})r. However, for z^ G H'^{^) n H^jj^n.), the trace theorem asserts 
that a^{z^) £ L (F) and (6, {(Tiy(z;^)})r coincides with its extension to a duality pairing. From 



the weak form of the Volterra dislocation problem ( 11 ) it moreover follows that ar(u, z^p) = l{zip) 



The Galerkin-orthogonality property of the WSM approximation (23) then yields the following 
identities: 



[u - M/i, ip)n = l{z,^ ~ Vh) ~ ariuh, z^ - Vh) + (6, {a^{vh)})r - (&, {a,y{z^)})T 
= ariu - Uh, z^ - Vh) - {b, {a^{z^)} - {<T„{vh)})r 



(58) 



which are valid for all Vh G Vh- 

For the first term in the final expression in ( 58 ) we can construct an appropriate bound without 
further digression. However, the second term, pertaining to the difference between the average 
traction oi z <E H (fl) n Hq x)(r2) and the average of the direct evaluation of the traction of the 
finite-element function Vh S V/i, is more difficult to estimate. The essential complication in con- 
structing an esti mate , is that the direct evaluation of the traction of the Galerkin finite-element 
approximation of (56) in Vh does not coincide with the weak formulation, because Vh ^ -H^divcrC^)) 
see also [2U]. Approximation results are available for the weak formulation of the traction (see, for 
instance, [5) but, to our knowledge, not for the direct formulation. 

To estimate the second term in ( 58 ) , we will use an auxiliary Nitsche-type approximation |15j to 



the dual problem. We consider the broken approximation spaces Vh according to (19). Moreover, 
we define V{h) := {H^{n \ F) n Hljj{n \ F)) -I- Vh- We equip V{h) with the mesh-dependent 
inner product 

iu,v)^^^-^^ar{u,v) + CY.h-HM,[[v]]), (59) 

sESh 
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and the corresponding norm || • Hv^c^-)- We note that the embedding of H'^{fl \ T) into V{h) 
is continuous, i.e., |1w||-(>//j1 < '^^||w||2,a\r for all v S i?^(Sl\r). We define the bilinear operator 
ar : V{h) x Vh — > K according to: 



ar(u, v) = ar{u, v) - ([[^;]], {cT,(^.)})^ + 9 ^ ([[u]], [[vl 

seSh 



(60) 



with 9 a suitable constant. For a suitable choice of the constants ( and 9, the bilinear form in the 
left member of (61 1 is continuous and coercive on Vh x V^ and ar{u, •) : Vh ^^ U^ represents a 



continuous linear functional for all u e V{h). Let z^^ e V/i denote the solution to the following 
Nitsche-type projection problem: 



z^h ^ ^h ■ ariz^h^'^h) = ariz^,Wh) Vw^ G Vh 



(61) 



By virtue of the continuity and coercivity of ar on Vh x Vh and the continuity of ariz^p, •) on Vh, 
the projection problem (61 1 defines a unique element z^^^ e Vh which satisfies 



^Vh\lV{h) 



< '€ inf 



2tp 'Wftllv(h) 



< 



inf |z;p-w/i|i,n 
wheVh 



(62) 



The first inequality in ( 62 ) is a straightforward consequence of Cea's lemma. The second inequality 



follows from V h C V h-, the continuity of functions in V h and (14a I 



By adding a suitable partition of zero to the second term in ( 58 ) and applying the triangle 
inequality, we obtain: 



|(5,{cr,.(z^)} - {<Jy(yhy\)v\ < \{b,{a^{zp ~ Zp^)})r\ + \{b,{a^{vh - Zpf^)})i 



(63) 



For the first term in the right- member of (63), we derive from (61 1 and (62): 

\{b,{(T^{zp -z^^)})r| = {b- [[wh]],Wi^{z<p - z^i^)})^ + ar{wh, z^ - z^,^) 



seSh 



h-\[[whl[[z^h]])^ 



(64) 



< ||^-[[wft]]||r(IIK(^v)}llr + IIK(2^J}llr) 
+ CA\wh\i,n\r\z,p - Zpf^\i,a\r 

+ ^Uwh]]\\j-^^-p\\[[zphMTH,r 

for all Wh G Vh- The trace theorem implies ||{o'i/(2i^)}||r < '^\\z,p\\2,n < "^^Hiysllo- The continuity of 
Or and ar implies that ||{o'i/(zi^^)}|jr < '^^\\ZiPh\\v(h)- Moreover, by virtue of (62) and the continuity 
of the embedding oiH'^{n) into V{h) it holds that ||2:^,Jlv.(,,) < '^IKUi^h) ^ "^ll^^lU^n < "^Mn- 
Hence, we have \\{a,{zp)}\\r + ||K(2^J}||r < '^M\n- Noting that Th = 7r{Vh) D Hl^^iV) 
coincides with {[[w;;]] : Wh G Vh}, for all A^ G Th there exists a Wh G Vh such that [[ui/i]] = Xh 
and |w/i|i,n\r < "^11^111/2, r- From (62)-(64) we deduce: 



\{b,{ai,{z^ - z^fj})r\ < '^[\\ip\\n\\b- \h\\j. + 



(||A/>||i/2,r + ||A/i||r;,,r) inf \zp ~ Wh\i,n\r 

WhEVh 



20 



for all A;, e Tfi, with || • WthX according to (43). For the second term in the right- member of (63 1 



we deduce from (47), the triangle inequality and (62): 

1(5, {a,K - 2^J})r| < 2-i/2c^CrMi/2||6|!r„rbh - ^^Ji,n\r 

<''^\\b\\rh.ri\z,p-Vh\i,n+ inf \z^ - Wh\i,n) 



(65) 



with ^ independent of h. Collecting the results in (58)-(65l and taking the infimum with respect 
to Vfi and A^, one obtains the estimate in (54). D 



Under slightly stronger conditions on the regularity of the slip distribution 6, we can derive from 
Theorem [4] a straightforward characterization of the asymptotic approximation properties of WSM 
in the L -norm for quasi-uniform meshes in the limit as ft. — >■ 0: 

Corollary 5. Assume that the conditions of Cor ollary [3| and Theorem^hold and, moreover, b G 

1 1/2 I, I 

H (r)ni?Q (r). Let u denote the solution to (11 ) and let Uh G V h denote the WSM approximation 
according to (23). It holds that 

\\u~UH\\n<'^h^'^ 

as h ^)- 0. 



(66) 



Proof. Based on the estimate ||2:;p||2,n < '^|i<(5||o, it follows from standard interpolation theory in 



Sobolev spaces (see, for instance, [HE]) that infi,,^gVh \z 
H\r)nHl/^{T),wehavemix,GT,\\b- XhWr <'<^h\\b\\,T- 
holds that |lA;;||i/2,r < "^ and ||A;!; 



Vh\i,n < '^h\\ip\\Q. Similarly, for b in 



turn, it follows from Theorem^ that 
into (54 1, we derive: 



u 



Setting A^ 



arg inf 



AhSTfe 



||&-A^||r, it 



< '(^h-^/'^. From \b2\ we obtain \\b\\rnX < ^^"^^^ and, in 

0. Inserting the above estimates 



\u-Uh\\i-i= sup 

¥)ei^(o)\{o} 

< sup 

vei^(o)\{o} 



w/Ji,o <'^/i~i/2 as ft 

\{u~Uh,Lp)n\ 

\W\\n 

-^(^h-^/''h\Mn + ^hMn 



1 + h 



-1/2 



)hM 



as ft — )• 0. The estimate in (66 1 then follows straightforwardly by combining terms. 



(67) 
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Remark. The reinforced regularity condition b E H (F) n Hq (F) on the slip distribution ensures 
that infAjjGTh \\b ~ Xh\\r < '^h as ft — >■ 0. Noting that the weaker estimate infA^gTh \\b — Xh\\r < 
■^ft^/^ suffices to obtain the result in Corollary[5J one is lead to question whether the reinforced reg- 
ularity condition on b is actually necessary. If the condition is dismissed, however, an interpolation 
estimate in the fractional Sobolev space Hq (F) is required. Interpolation estimates in fractional 
Sobolev spaces are technical (see, for instance, [S]) and the particular result required here is to our 
knowledge not available. 

It is noteworthy that the asymptotic convergence behavior according to Corollary [5] is consis- 
tent with the notion that the continuous WSM approximation incurs an 0(1) error, pointwise, in 
the 0(ft) neighborhood of the discontinuity composed of the intersected elements as ft — ?> 0. In 
particular, denoting by Th the union of the elements for which the intersection of the fault with the 
closure of the element is non-empty, it holds that 

||l||r, = (measw(F;,))'/' = 0(fti/2) (68) 

as ft — ?> 0, with measAr(F/i) the A^-Lebesgue measure of T^. 
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4.3 Local approximation properties 

Next, we consider the local approximation properties of WSM, i.e., on f2 excluding a small neigh- 
borhood of the fault. To this end, we will relate the WSM approximation to the standard Galerkin 
approximation associated with a locally supported lift, outside the support of the lift. 

For arbitrary e > 0, let k^ :— {x ^ H. : dist(a;, k) < e} denote the open e-neighborhood of the 
dislocation fault. We consider a lift 4 G Hljj{fl\T) of b with compact support in k'^ . Let Uh & Vh 
denote the Galerkin approximation of the continuous complement of the solution with respect to 



the jump lift ib according to ( 17 1. As a straightforward consequence of Cea's lemma, it follows that 



u/i + lb is endowed with the quasi-optimal approximation property: 

\\u-{uh+ib)\\i.n\r<'^ inf \\u - {vh + £b)\\i.n\r (69) 

A meaningful characterization of the local approximation properties of WSM is therefore provided 
by an estimate for the deviation \\uh + (-b~ Uh\\i,Q,\>c'' between the WSM approximation Uh and the 
local-lift-based approximation Uh + £b- Note that >f'^ D supp(^f,) implies that £b vanishes on il \ >f ^ 
and, hence, the estimate pertains to the deviation between Uh G Vh and Uh G Vh- 

The characterization of the local approximation properties of WSM in the iJ^-norm in The- 
orem |6] below is based on the interpolation of a particular extension of functions in Vh onto fi^ . 
Specifically, we consider an extension corresponding to the operator E : Vh -^ H^{fl): 

E{vh) ^Vh in rj \ x^ (70a) 

-diY a{E{vh)) =0 in ^^ (70b) 

in the sense of distributions, and its optimal approximation in the norm defined by ar(-, ■): 

Eh{vh) = arg min ar{E{vh) - Wh, E{vh) - Wh) (71) 

where V/j^j^e := {vh € Vh '■ supp(D;i) C x^} ^ denotes the class of approximation functions of 
which the support is confined to x^ . Note that Vh,x'^ 7^ implies e > 'rfh. 

Equations (70a[) and (71) imply that E{vh) and Eh{vh) coincide with Vh outside x^ . Inside x^. 



the extension _E(u;i) is defined by the homogeneous elasticity problem (70b I. Fromi5(w;i) g H (J7) it 



follows that E{vh) is continuous at dx^ in the trace sense, which implies that (70b) is complemented 



by the Dirichlet boundary condition E{vh) = Vh on dx^ . By virtue of the smoothness condition on 



the elasticity tensor in (27), the extension operator according to (70) exhibits an interior regularity 
property on x^ . In particular, for all Vh € Vh it holds that E{vh) G H^^^i^x'^) (i.e., (j)E{vh) € 
H (>f^) for any (j) G C°°(51) with compact support in x^) and the following estimate holds for each 
open subset K (s^ x'^: 

\\E{vh)\\2,K < '^\\E{vh)\\^^ -ivh e Vh\ (72) 



see, for instance, [51 §6.3.1] or [T31 §2.3] for further details. It is important to note that estimate ( [72| 
in conjunction with the trace theorem implies that 

\\{a,[E{vh)))\\^<^E{vh)\\^^^^ yvh(.Vh, (73) 

for some ^ > independent of w/i, provided that there exists an open subset K (S: x^ such 
that X d K. This provision implies that the dislocation must be properly contained in f2. The 
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proof of Theorem l6] involves an estimate of the difference between the average traction of E{vh) and 
its approximation Eh(vh) according to (71 1. In general, we can estimate (6, {(Tv[E[vh) — Eh{vh)y\)v 
in the same manner as in the proof of Theorem Bj The bound ( 73 1 in combination with the optimal- 
approximation property of Eh{vh) in (71 ) however suggests that in this case a sharper estimate can 



be established. The derivation of such a refined estimate is intricate, however, as (&, {o',y(-)})r is 
unbounded on H^{^) and the estimate involves the difference between E{vh) and Eh{vh), which 
reside in different subspaces of ff ^(fi). We therefore formulate the refined estimate in the form of 
a conjecture: 

Conjecture 1. Assume that there exist an e* > Q and an open subset K <! x'^ such that x (Z K . 



Consider a sequence of approximation spaces V-h, the extension operator according to (70) and its 
approximation according to \71^ . For all e > e* and all Vh S Vh, it holds that 

\\{a,iEivh) - Eh{vh))}\\^ < ^Eivh) - Eh{vh)\\i,^e (74) 

for some constant "^ > independent of h. 

Theorem l6] below presents a general characterization of the local approximation properties of 
WSM, independent of Conjecture [ll and a refinement, which is contingent on the conjecture. 

Theorem 6 (Local approximation properties of WSM in the H -norm). Assume that there exist 
an e* > and an open subset K (s x^ such that x C K . For arbitrary e > e* , let £b € Hq x)(f^ \ T) 
denote a lift of b such that supp(^b) C x"^ . Given a sequence of approximation spaces V-h, for each 
h G H let Ufi denote the approximation of the continuous complement corresponding to £b in [P, 



and let Uh G V t denote the WSM approximation according to (23). Consider the extension operator 
in (70) and its approximation in Vh according to (71). It holds that 



\\uh +4- Wh.||?,fi\^e < '^UbWi^x^xrWEiuh -Uh) - Eh{uh - Uh)\\i^^e 
+ ^\{b,{ay[E{uh -Uh) - Eh{uh - M/i))})p| 
For some constant ^ > independent of h and e. If, in addition, ConjectureUl holds, then 

II- , fj II ^ c^rw/j w , iij,ii \ ■ ( \\E{vh) - Wh\\i.>t^ 
\\Uh+ib-Uh\\i^n\>,- <^[\\h\\i.>c-\r + \\b\\r) sup mf — ■ — 

for some constant "^ > independent of h and e. 



\vh\\i,n\> 



(75) 



(76) 



Proof. We use the condensed notation Eh := Eh{uh — Uh). Noting that (71) implies that Eh 
coincides with Uh — Uh on fl\ x^ and that supp(i!f,) C >f^, it holds that 



\Uh 



h - uh\\x.a\>c'- = \\uh - uh\\\.n\> 



< WE, 



h\ l.fl 



(77) 



Using the Poincare inequality ( 15 ) and the strong positivity of ar(-, •) according to ( 14b I, we obtain 
the following bound: 

\\Eh\\ln<i^ + Cp)cA'\ar{Eh,Eh)\ (78) 



By introducing a suitable partition of zero, we derive from the WSM formulation (23) and the 



lift-based Galerkin approximation (17) 



ar{Eh,Eh) = ar(u/i - Uh,Eh) + ar{Eh - (uh - Uh),Eh) 

= (6, {crj,(£'/i)})p -ar(4,^h) + ar{Eh - {uh -Uh),Eh) 



(79) 
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By virtue of ( |70b ), ar{vh,E{uii — Uh)) vanishes for all Vk £ Vh,x^- The optimality condition 
associated with (71 1 therefore reduces to: 



arivh,Eh) = 



yvh e Vt 



(80) 



Noting that Ef^ — (u^ — Uh)& 'Vh,3<^ : Equation ( 80 ) implies that the right-most term in the ultimate 



expression in ( 79 1 vanishes. By virtue of the interior regularity of E{uh — u^), the following identity 
holds: 

aT{eb,E{uh~Uh)) - {b,{a,{E{uh-Uh))})^ (81) 



Let us note that (81) corresponds to the (admissible) identification of the dualit y pa ir ing between 

b e ify^(>f) and {a^{E{uh - w/i))} G H^^^^{>c) to an L^ inner product. From (79H(81) and the 
triangle inequality we then deduce that: 



\ar{Eh,Eh)\ < \ar{eb,E{uh - Uh) - Eh)\ + \{b,{a^{E{uh ~ u,,) ~ Eh)})i 



(82) 



We recall that E{uh ~ Uh) — Eh vanishes on 51 \ >f^. Estimate (75) then follows straightforwardly 
from (77 1, (78 1, (82) and the continuity of ar(-, •) according to ( 14a[ ). 

To prove the auxiliary assertion (76), we note that \\E{vh) — Eh{vh)\\i,Q in fact depends only on 



the trace of u^ on dx'^ and, by the trace theorem, it holds that \\E{vh) — Eh{vh)\\i,n < '^||w?ilii,o\> 
Therefore, if Conjecture fl] holds, we obtain from (JTSl): 



\\uh -'"/i||?,n\>c= < '^(Il4||i,>c=\r + ll&l|r)||£^(u/i - u/i) ~ Eh{uh - Uh)\\^^ 

\E{vh) - Eh{vh)\\ 



<'^(ll4||i..=\r 



l^llr) sup 



(83) 



\Uh ~Uh\\i,n\> 



Estimate ( 76 1 follows directly from the identity in ( 77 ) , the ultimate bound in ( 83 1 and the definition 
of Eh in ^. 
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Theorem p^ essentially implies that if the lift-based approximation (17) displays optimal global 
convergence, then the WSM approximation (23) displays optimal local convergence. 



Corollary 7. Assume that Conjectured^ holds and that there exist an e* > and an open subset 
K (s^ >if such that h C K. Assume that Aijki G C^^^{^) for some integer k > and that fl is 
convex or of class C^'^'^ . Assume that b admits a sufficiently smooth local lifting, in particular, that 
there exists an ib such that {cr^(ib)} — 0, div cr(£b) S H (57) and supp(^b) C >if . Let V-u denote 
a sequence of H (il)- conforming piecewise polynomial approximation spaces of degree p > 1 with 



approximation property (18). Let u denote the solution to the Volterra dislocation problem (11) and 



let uu denote the sequence of WSM approximations { 23 ) corresponding to the approximation spaces 
V-H- It holds that 

\\u-uh\\i,n\>.^<'^h' (84) 

OS /i — >■ with ^ > independent of h and I = min{p, fc + 1}. 



Proof. Note that the following integration- by-parts identity holds for all v G H'^{n): 



ar{ib,v)= / w-(T„(4)- / w-divcr(4)= / v ■ {a^{eb)} 

Jd{n\r) Jn\r Jr 



V ■ div cr( 



(85) 



o\r 
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The second identity follows by rearranging terms and applying (|8|). The first term in the ultimate 



expression in (85) vanishes by virtue of the conditions on It,. Moreover, because l{v) = (/, w) — 0, 



the right member of (12 1 corresponds to 



l(v) - ar%,v) = - I wdivcr(4) (86) 

Jn\T 

By virtue of div cr(£b) € H [Vl) and the conditions on the elasticity tensor and the domain, it 



then holds that the continuous complement u in (12) resides in H ^^(rj); see, for instance, [5| 
Theorem 6.3.5]. The quasi-optimality of the Galerkin approximation (17) in combination with 
the optimal approximation properties of Vh according to (18) implies that ||u — u/i||i.t2 < "^hK 



Moreover, the conditions on the elasticity tensor imply that E{vh) G H^^'^{>€^) and, in turn, (18) 



yields: 

sup inf \\EM_p^!4l^^^j^i (87) 



The estimate ( 84 ) then follows from 

\\u ~ Uh\\i,n\^e < ||m-u/i|1i,o + \\uh + ib - Uh\\i,n\>,^ 
and Theorem |6l D 

5 Numerical results 

We will assess the approximation properties of the Weakly-enforced Slip Method on the basis of 
three different test cases. All three test cases are designed to have analytical results available. This 
allows us to compare the exact solution of boundary value problem (Is]) with the WSM solution 



( 23 ) and study the behavior of errors and convergence under refinement of the finite finite-element 
mesh. 

The considered test cases are: 

I. a two-dimensional infinite domain loaded in plane strain, with a straight, finite dislocation 
and smooth slip distribution, 

II. a three-dimensional semi-infinite domain with traction-free surface, a planar, non-rupturing 
dislocation and constant slip, and 

III. a three-dimensional semi-infinite domain with traction-free surface, a planar, surface rupturing 
dislocation and constant slip. 

All three test cases are in principle set on infinite domains. To make the problems amenable 
to treatment by the finite-element method, we truncate the domains, and restrict the analyses to 
suitably large, but finite, neighbourhoods of the dislocation. In order to focus on the treatment 
of the dislocations, boundary truncation errors are controlled by constraining the finite-element 
approximation to the analytical solution at the artificial lateral boundaries. 

Let us note that of the three test cases, only the first one is by design in full accordance with 
the theory developed in Section |4] Test case II and III feature a piecewise constant slip to match 
available analytical solutions, violating the condition of Lemma [T] which states that the slip can be 
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V 



Figure 4: Schematic overview of test case I. The dashed hne marks the computational domain, which is 
a truncation of the actual (infinite) domain. A unit-length dislocation x is placed at an arctan(3/4) angle 
with the horizontal axis at the center of the computational domain. 



lifted into HIt^{Vl\T 



i.e., h € -ffp (-^)- Test case III moreover considers a rupturing fault and, 
accordingly, it has nonzero slip at the intersection with the domain boundary. The results below 
convey that the main results of the theory nonetheless uphold, suggesting that the conditions under 
which they apply can be relaxed. 



5.1 Test case I: 2D plane strain 

The first test case is a line dislocation in plane strain in an infinite two-dimensional isotropic domain. 
The dislocation is straight and of unit length, with a smoothly varying slip that is tangent to the 
dislocation line, making it a pure shear dislocation. Figure [4] shows a schematic of the computational 
setup. 

Introducing an arclength coordinate ^ along the fault, we consider a smooth, piecewise quadratic 
slip distribution according to: 




2 ^ "^ ^ fi 



otherwise 



(89) 



The dislocation, which corresponds to the support of the slip, is located in the interval ^ € [— |, |]. 
The slip is symmetric, with zero displacement at the tips and smoothly opening and closing. We set 
the scaling factor feo = 0.1. The value of bg is however non-essential on account of the linearity of the 
problem. An analytical solution to this problem can be constructed for an infinite, homogeneous, 
isotropic domain, by adapting the well known solution for edge dislocations [10] . Introducing a 
perpendicular coordinate 77, the resulting displacement along the (^, 77) coordinates is expressed in 
terms of Lame parameters A and /i as 



M(f , 77) = 60 [6C/(e - i 77) - 18t/(e - i, 77) + 18C/(e + i, 77) - 6C/(e + i, 77)] 



(90) 



where 



Ui^.il) 



M C2 _ 2A+A1 2 



logv/^^T^ 



A+2m 



A+2/j' 
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Figure 5: Side by side view of tlie discontinuous exact solution of test case I (left) with the continuous 
WSM approximation on a 16 x 16 mesh and linear shape functions (right). Color represents displacement 
magnitude on a log scale. Mesh lines on the left are introduced to indicate displacements and to facilitate 
comparison. 



A finite sized problem is obtained by truncating the domain to J7 = (—1, 1)^ and introducing a 
Dirichlet condition at the boundary dQ ~ V with data corresponding to the exact solution. The 
fault is located in the center of the domain at an angle of arctan(3/4). All numerical results are 
generated for Lame parameters A = 1, ^ = 1. 

Figure [5] shows the exact solution of Eqn. (901 (left) in the deformed configuration, side by side 
with the WSM approximation for linear shape functions on a 16 x 16 mesh (right). The colors 
indicate displacement magnitude on a log scale. The grid-line pattern represents the distortion of 
a mesh that is uniform in the undeformed configuration. The lines coincide with element edges for 
the finite element computation on the right, and are matched on the left to facilitate visual com- 
parison. Because the test case is symmetric, we expect the images to be approximately rotationally 
symmetric. We observe that the rotational symmetry breaks at the dislocation, where the exact 
solution is discontinuous, while the WSM approximation is continuous throughout the domain. 

The global symmetry of the displacement pattern of Figure [5] indicates that the error induced by 
the WSM approximation diminish with distance from the dislocation. To quantify this observation. 
Figure [6] shows the errors for a 16 x 16 mesh (left), and for a sequence of meshes along the line 
B-B' perpendicular to the fault (right). The latter shows that for every mesh the error decays 
exponentially with increasing distance from the fault, with the exception of a narrow zone in the 
vicinity of the dislocation, where the error displays a transition to the constant error ^bo at the 
dislocation. One may observe that the error at the dislocation is independent of element size. 
Away from the dislocation, errors are seen to decrease with element size, approximately reducing 
by 10*^'^ « 4 whenever the mesh width is halved. This error reduction provides a first indication 
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Figure 6: Absolute value of the displacement error for the 16 x 16 mesh result of Figure pi (left), side by 
side with error evaluations along the line B-B" for a sequence of meshes of increasing density. The errors 
shown on the left thus correspond with the 3rd (red) curve on the right. 



that local L^ convergence is optimal at 0(hP'^^) as ft. — )• 0. 

Convergence of the WSM approximation under mesh refinement is further examined in Figure [Tj 
which displays the L -norm and H -norm of the error for both linear (p = 1) and quadratic {jp = 2) 
shape functions. The curve marked 'global' shows the error integrated over the entire computational 
domain. The curve marked 'local' shows the error integrated over the domain excluding an 0.1- 
neighborhood of the dislocation, corresponding to the dotted area in Figure [6J We observe that for 
both linear and quadratic approximations, the global H -norm of the error diverges as 0{h~^'^) 
as /i — >■ 0, while the global L -norm of the error converges as 0{h~^^^), independent of the order 
of approximation. This asymptotic behavior is in accordance with the estimates in Corollaries |3] 
and 5l Figure M moreover corroborates that the local iif ^-norm of the error converges as 0(hP) as 
/i — >■ 0, in agreement with the estimate in Corollary [t] The local i^-norm of the error, for which 
no theory was developed, also displays an optimal convergence rate of 0(/iP~'"^) as /i — ?► 0. 



5.2 Test case II: 3D traction-free halfspace 

The second test case is a planar dislocation buried in a semi-infinite, three dimensional, homoge- 
neous, isotropic domain with a flat, traction-free surface. The dislocation is a rectangular plane, 
the sides of which are displaced over a constant distance in both strike and in dip direction, such 
that in geodetic terms the setting is that of an oblique left-lateral thrust fault. Figure [8] shows a 
schematic of the computational setup. Analytical solutions to this problem have been derived by 
Okada [T6l|T7j. We use homogeneous Lame constants A = 1, /.t = 1 for all subsequent computations. 
Computations are performed on a truncated domain 51 = [— 1, 1] x [— 1, 1] x [— 1, 0], with Dirichlet 
conditions enforcing the Okada solution at the five truncation planes. Figure [9] shows displacements 
along the intersection planes A-A' and B-B' indicated in Figure [8] The topmost figure shows the 
exact solution according to Okada's equations. The middle and bottom figures display the WSM 
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Figure 7: Mesh convergence of WSM applied to test case I, siiowing tlie L^-norm (top) and //^-norm 
(bottom) of tiie error for linear (left) and quadratic (right) shape functions. The markers corresponds to 
mesh sequences of {4 x 4, 8 x 8, ... , 128 x 128} elements are considered. The global error is computed by 
integration over the entire computational domain Q, the local error by integration over the >0.1 distance 
exterior around xr, bounded by the dotted line in Figure [HI 




-1,1] X [-1,1])^^-/ 1,0] 



> X 



Figure 8: Schematic overview of test case II. A 1 x 3~^" sized dislocation plane is positioned at 15° 
strike and 30° dip, spanning a 0.25-0.75 depth range. The surface is traction free. Dashed faces mark the 
computational domain, which is a truncation of the actual (semi-infinite) domain. The fault is represented 
by the B-B' plane, where at x the medium is dislocated by a constant 0.2 displacement in strike direction and 
a constant 0.1 displacement in dip direction. The perpendicular intersection plane A-A' serves visualization 
purposes only; see Figures [91 and 12 
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Figure 9: Cross-section views of the displacement field of test case II, corresponding to the A-A' and B-B' 
intersection planes of Figurels] exact displacement as derived by Okada [17] (top) and WSM approximation 
on a 16 X 16 X 8 mesh with linear (middle) and quadratic shape functions (bottom). Colors indicate the 
displacement magnitude on a logarithmic scale analogous to Figure [5| Arrows indicate direction only. 
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Figure 10: Mesh convergence of WSM for test case II. The figures are analogous to Figure ItI with the 
addition of the L^-norm of the error in the displacement at the traction-free surface. 



approximation of the displacement on a 16 x 16 x 8 mesh, for hnear and quadratic shape func- 
tions, respectively. The displacement fields in the WSM approximations are seen to be continuous 
everywhere, though they become highly irregular at the dislocation. At further distances from 
the dislocation, however, the approximations exhibit very good agreement with the exact solution, 
despite the coarseness of the considered mesh. 

Figure [TO] examines convergence of the global and local norms of the error under mesh refine- 
ment. The results confirm the OihT^I"^) divergence and the 0(/i^/^) convergence of the global 
H -norm and global h -norm, respectively, in agreement with the estimates in Corollaries^ and pi 
FurthermorCjWe observe 0[W) convergence of the local iJ^-norm in accordance with the estimate 
in Corollary M and 0(hP'^^) convergence for the local L -norm. Figure 10 moreover displays the 
h -norm of the displacement error at the surface, which exhibits an optimal convergence rate of 
0(hP^^). It is to be noted that the agreement of the observed convergence rates for the global 
H -norm and h -norm and the local H -norm with the estimates in Section W^ is non-obvious, as 
a piecewise constant slip cannot be lifted into if g p(il \ T), and the test case under consideration 
therefore fails to meet the conditions imposed in Section |4J 

5.3 Test case III: 3D traction-free rupturing halfspace 

The third test case is in everything equal to the second, except that the dislocation is now extended 
in vertical direction to span a depth range of 0-0.75, and the dip slip direction is reversed to form 
the geodetic equivalent of a rupturing reverse fault; see Figure [TT] Figure [12] shows the cross 
section displacements, and Figure [13] the norms of the displacement-error under mesh refinement. 
We observe that the global (resp. local) H -norm and h -norm of the error are again proportional 
to hr^l"^ and h^l'^ (resp. hP and hP^^)^ respectively. In addition to the L -norm of the error in 
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Figure 11: Perspective projection of a WSM computation of test case III, using linear shape functions on 
a 128 X 128 x 64 mesh. Deformations are amphfied by a factor 2 for visuahzation purposes; colors represent 
displacement magnitude on a log scale. The orientation aligns roughly with that of Figure [8] A spherical 
cut-out exposes part of the interior of the domain, revealing a zone of large displacements local to the 
dislocation. 



the displacement field at the traction-free boundary, Figure 13 also presents the local L -norm of 



the error at the traction-free boundary, i.e., the error on the surface excluding an 0.1-neighborhood 
of the dislocation. The global L -norm of the surface error converges with an asymptotic rate 
of 0{h^^^). This suboptimal convergence behavior is caused by the fact that in this case the 
discontinuity reaches the surface. The local L^-norm of the error at the surface again display 
optimal convergence at a rate of 0{hP'^^). It is to be noted that the convergence results for the 
global H -norm and L -norm and the local H -norm agree with the estimates in Section |4J despite 
the fact that the analysis in Section [4] is restricted to non-rupturing faults. 

6 Conclusions 

To solve Volterra's dislocation problem by standard finite-elements techniques, the dislocation is 
required to coincide with element edges. This requirement links the finite-element mesh with the 
fault geometry, which prohibits the reuse of computational components in situations where multiple 
geometries have to be considered. In particular, it renders the finite-element method infeasible in 
nonlinear inversion problems. 

To overcome the problems of standard finite-element techniques in nonlinear inversion processes, 
in this paper we introduced the Weakly- enforced Slip Method (WSM), a new finite-element approx- 
imation for Volterra's dislocation problem in which the slip discontinuity is weakly imposed in the 
right-hand-side load functional. Accordingly, the bilinear form in the formulation and, hence, the 
stiffness matrix are independent of the fault geometry. The method is summarized by the following 
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Figure 13: Mesh convergence of WSM applied to test case III. Tiie figures are analogous to Figure fO 
except that the i^-norm of the error on the surface is now indicated by 'global surface' and additionally a 
'local surface' L^-norm of the error has been plotted, excluding an 0.1-neighborhood of the dislocation. 



weak formulation: 



u^Vh'- a(u, v) = - b ■ {<7„{v)} yv eVh' 



The stiffness matrix depends on properties of the continuous domain only, namely, variations in the 
elasticity and topology and geometry of the domain, and remains independent of fault geometry. 
Fault dependence manifests in the right-hand-side load functional only. The load functional is 
formed by integrating over the fault, which is allowed to cut through elements. The integration 
along the fault is a non-standard operation in finite-element methods, but we expect that it can 
be incorporated in most existing finite-element toolkits with minor effort. We further note that 
no approximations are required regarding fault geometry and slip distribution, unlike lift-based 
methods, which require a parametrization for both. 

We established that as a consequence of the continuous approximation in WSM, the dislocation 
is not resolved. Accordingly, the WSM approximation displays suboptimal convergence in the L - 
norm under mesh refinement. In particular, the L -norm of the error decays only as 0{h}'^) as the 
mesh width h tends to 0, independent of the order of approximation. Furthermore, the H -norm 
of the error generally diverges as 0{h~^'^) as /i — > 0, independent of the order of approximation. 
In addition, we however proved that WSM has outstanding local approximation properties, and 
that the method generally displays optimal convergence in the H -norm on the domain excluding 
an arbitrarily small neighborhood of the dislocation. In particular, for any e-neighborhood >f^ of 



the dislocation, the 



1,S1\> 



norm of the error in the WSM approximation generally converges as 



0{hP) as ft. — >■ 0, with p the polynomial degree of the finite-element space. 

Numerical experiments in 2D and 3D were conducted to verify and scrutinize the approximation 
properties of WSM. The asymptotic error estimates for the global i?^-norm, the global i^-norm 
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and the local Jf ^-norni of the error were confirmed in all cases, despite the fact that two of the test 
cases violate some of the conditions underlying the asymptotic error estimates. In particular, the 
numerical results indicate that the error estimates extend to rupturing faults, where the dislocation 
fissures the traction-free surface. The numerical experiments moreover conveyed that the WSM 
approximation displays optimal local convergence in the L -norm, i.e., the || • ||{2w=-norm of the 
error decays as 0{h^^'^) as ft. — )■ 0. For the non-rupturing-fault test case, the error in the surface 
displacement converges optimally in the X^-norm at 0(/iP+^). For the rupturing- fault test case, the 
L -norm of the error at the surface converges at 0(/i"'^/^), while the local L -norm excluding a small 
neighborhood of the dislocation again converges optimally at 0{N''^^). Overall, the approximation 
obtained via WSM is very well behaved, and the method proves remarkably robust. 

Given the compelling properties of WSM one might be tempted to seek application in other than 
our intended field of tectonophysics. Dislocation plasticity comes to mind as one heavily relying on 
superposition of elastic dislocations. For many problems, however, the location of the fault will be 
known a-priori, in which case there is no reason to avoid strong imposition. Moreover, often the 
internal stress is an important quantity to be resolved, which with WSM suffers from inaccuracies 
in regions close to the fault. In tectonophysics the primary observable is the displacement of the 
free surface, which WSM is very well capable of resolving. 

We believe that WSM can play an important, if specific, role in the application of tectonic fault 
plane inversions. Being able to precompute the stiffness matrix and a quality preconditioner, for any 
fault geometry and slip distribution under study, it remains only to integrate over the 2D manifold 
and solve the system. This makes it feasible to use finite elements in a direct nonlinear inversion. 
In the hands of geophysicists, this tool will allow all available in-situ knowledge to be made part of 
the forward model. We hope this will help to improve the accuracy of future co-seismic analyses. 
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